#Glasser region and label names
glasser.regions <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/HCPMMP_glasseratlas/glasser360_regionlist.csv")
glasser.frontal <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/HCPMMP_glasseratlas/glasser360_regionlist_frontallobe.csv")
glasser.snr.exclude <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/SNR/glasser_SNR_exclusion.csv")

depth.list <- c("depth_1", "depth_2", "depth_3", "depth_4", "depth_5", "depth_6", "depth_7")
ordered_depths <- c("depth_7", "depth_6", "depth_5", "depth_4", "depth_3", "depth_2", "depth_1")
depth_colorbar <-  c("#004A38FF", "#006B63FF", "#008FA7FF", "#6992CC", "#A29DC4", "#C0BFE3", "#E2CFE5")
depth_colorbar_reverse <-  c("#E2CFE5", "#C0BFE3", "#A29DC4", "#6992CC", "#008FA7FF", "#006B63FF", "#004A38FF")

setwd("/Volumes/Hera/Projects/corticalmyelin_development/gam_outputs/developmental_effects/") #output from /gam_models/fitgams_glasserregions_bydepth.R

files <- list.files(getwd(), pattern = "age.RDS") 

#read in files and assign to variables
for(i in 1:length(files)){
  
  Rfilename <- gsub(".RDS", "", files[i])
  
  x <- readRDS(files[i]) 
  assign(Rfilename, x) 
  }

gam.results.alldepths <- list(depth1_gamstatistics_age, depth2_gamstatistics_age, depth3_gamstatistics_age, depth4_gamstatistics_age, depth5_gamstatistics_age, depth6_gamstatistics_age, depth7_gamstatistics_age)
names(gam.results.alldepths) <- list("depth_1", "depth_2", "depth_3", "depth_4", "depth_5", "depth_6", "depth_7")

gam.smoothestimates.alldepths <- lapply(gam.results.alldepths, '[[', "gam.smoothestimates.df" )
gam.smoothestimates.alldepths <- lapply(gam.smoothestimates.alldepths, function(depth){
  depth <-  depth[!(depth$orig_parcelname %in% glasser.snr.exclude$orig_parcelname),] #exclude low SNR parcels
})

gam.smoothestimates.alldepths.long <- do.call(rbind, gam.smoothestimates.alldepths)
gam.smoothestimates.alldepths.long <- gam.smoothestimates.alldepths.long %>% mutate(depth = as.factor(substr(row.names(gam.smoothestimates.alldepths.long), 1, 7))) #add depth as a column
gam.smoothestimates.alldepths.long$depth <- factor(gam.smoothestimates.alldepths.long$depth, ordered = T, levels = ordered_depths)

#Function to generate a ggseg cortical plot to indicate the location of region of interest 
plot.cortical.region <- function(region, hemi, side, save = TRUE) {
  
  regional.ggseg.plot <- glasser.regions %>% filter(orig_parcelname == region) %>%
  ggseg(.data = ., atlas = "glasser", mapping=aes(fill = orig_parcelname, colour=I("gray25"), size=I(.07)), position = c("stacked"), hemisphere = hemi, view = side) + 
  theme_void() + 
  labs(fill="") +
  scale_fill_manual(values = c("#395F94"), na.value = "white") +
  theme(legend.position = "none")
  
  regional.ggseg.plot
  if(save == TRUE){
  ggsave(filename = sprintf("/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure4_supplementary/%s_corticalmap.pdf", region), device = "pdf", dpi = 500, width = .85, height = .55)
  }
  return(regional.ggseg.plot)
}

#Function to plot age smooth functions in exemplar regions
plot.depth.smooths <- function(region, y_ticks, save = T){

  gam.smoothestimates.alldepths.long$depth <- factor(gam.smoothestimates.alldepths.long$depth, ordered = T, levels = depth.list)
  
  smooth.plot <- gam.smoothestimates.alldepths.long %>% filter(orig_parcelname == region) %>%
  ggplot(., aes(x = age, y = est, group = depth, color = depth)) +
  geom_line(linewidth = .4) +
  theme_classic() +
  xlab("\nAge (years)") +
  ylab("R1 Trajectory (zero-centered)\n") +
  scale_y_continuous(breaks = y_ticks) +
  scale_color_manual(values = depth_colorbar_reverse) +
  theme(
        legend.position = "none", 
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_text(size = 6, family = "Arial", color = c("black")),
        axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
        axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
        axis.line = element_line(linewidth = .2), 
        axis.ticks = element_line(linewidth = .2))
  
  smooth.plot
  if(save == TRUE){
  ggsave(filename = sprintf("/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure4_supplementary/%s_depthsmooths.pdf", region), dpi = 500, width = 2.5, height = 1.7)
  }

  return(smooth.plot)
}